Quantum computation and analysis of Wigner and Husimi functions: 
toward a quantum image treatment 
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We study the efficiency of quantum algorithms which aim at obtaining phase space distribution 
functions of quantum systems. Wigner and Husimi functions are considered. Different quantum 
algorithms are envisioned to build these functions, and compared with the classical computation. 
Different procedures to extract more efficiently information from the final wave function of these 
algorithms are studied, including coarse-grained measurements, amplitude amplification and mea- 
sure of wavelet-transformed wave function. The algorithms are analyzed and numerically tested on 
a complex quantum system showing different behavior depending on parameters, namely the kicked 
rotator. The results for the Wigner function show in particular that the use of the quantum wavelet 
transform gives a polynomial gain over classical computation. For the Husimi distribution, the gain 
is much larger than for the Wigner function, and is bigger with the help of amplitude amplification 
and wavelet transforms. We also apply the same set of techniques to the analysis of real images. 
The results show that the use of the quantum wavelet transform allows to lower dramatically the 
number of measurements needed, but at the cost of a large loss of information. 
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I.INTRODUCTION 

In the recent years, the study of quantum information 
has attracted more and more interest. In this field, 
quantum mechanics is used to treat and manipulate infor- 
mation. Important applications are quantum cryptogra- 
phy, quantum teleportation and quantum computation. 
The latter takes advantage of the laws of quantum me- 
chanics to perform computational tasks sometimes much 
faster than classical devices. A famous example is pro- 
vided by the problem of factoring large integers, useful 
for public-key cryptography, which can be solved with 
exponential efficiency by Shor's algorithm 2]. Another 
example is the search of an unstructured list, which was 
shown by Grover 3] to be quadratically faster on quan- 
tum devices. In parallel, investigations of the simulation 
of quantum systems on quantum computers showed that 
the evolution of a complex wave function can be simu- 
lated efficiently for an exponentially large Hilbert space 
with polynomial resources 4, 5, 6, 7, 8, 9]. Still, there are 
many open questions which remain unanswered. In par- 
ticular, it is not always clear how to perform an efficient 
extraction of information from such a complex quantum 
mechanical wave function once it has been evolved on 
a quantum computer. More generally, the same prob- 
lem appears for quantum algorithms manipulating large 
amount of classical data. 

In the present paper, we study different algorithmic 
processes which perform this task. We focus on the phase 
space distribution (Wigner and Husimi functions) jloL lnj 
These functions provide a two-dimensional picture of a 
one-dimensional wave function, and can be compared 
directly with classical phase space distributions. They 
have also been shown in llj, llJ to be stable with re- 



spect to various quantum computer error models. Dif- 
ferent phase space representation which can be imple- 
mented efficiently on a quantum computer will be ex- 
plored, first the discrete Wigner transform, for which an 
original algorithm will be presented, and then a Husimi- 
like transform, first introduced in this context in [l^ . 
Recent proposals Hi mile] gave methods to measure 
or construct Wigner and Husimi functions on a quantum 
computer, using for example phase space tomography. 
These method will be analyzed and compared with new 
strategies, in order to identify the most efficient algo- 
rithms. Different techniques will be tested in order to 
extract information, namely measure of an ancilla qubit, 
measurement of all qubits, coarse grained measurement, 
and the use of amplitude amplification p^. In addi- 
tion, we will analyze the use of the wavelet transform 
to compress information and minimize the number of 
measurements. Indeed, wavelet transforms 0, are 
used in a large number of applications involving clas- 
sical data treatment, in particular they allow to reach 
large compression rates for classical images in standards 
like MPEG. Quantum wavelet transforms have been built 
and implemented and it was shown that 

they can be applied on an exponentially large vector in 
a polynomial number of operations. Numerical compu- 
tations will enable us to quantify the efficiency of each 
method for a specific complex quantum system, namely 
the kicked rotator. In general, it will be shown that a 
polynomial gain can be reached with several strategies. 
Since a quantum phase space distribution can be con- 
sidered as an example of a two-dimensional picture, we 
discuss in a subsequent section the use of the same tech- 
niques to treat images encoded on the wave function of 
a quantum computer, in a way similar to what is done 




in classical image analysis. This for example could be 
applied to images transmitted through quantum imaging 



II. QUANTUM PHASE SPACE DISTRIBUTIONS 
FOR A CHAOTIC QUANTUM MAP 

Classical Hamiltonian mechanics is built in phase 
space, dynamics being governed by Hamilton's equation 
of motion. Classical motion can be described through the 
evolution of phase space (Liouville) distributions. On the 
other hand, phase space is a peculiar notion in quantum 
mechanics since p and q do not commute. A wave func- 
tion is naturally described in a Hilbert space, for example 
position alone or momentum alone. Nevertheless, it has 
been known since a long time that it is possible to define 
functions of p and q which can be thought as quantum 
phase space distributions. The most commonly used is 
the Wigner function 10], defined for the wave function 

of a continuous system by: 



W{p,q) 



Mq + ^rHQ-^)ck' (1) 



This function involves the two variables position q and 
momentum p in a symmetric way (although it is not 
immediately apparent in the formula Q), and shares 
some properties with classical phase space probability 
distributions. Indeed, it is a real function, and satisfies 
fW{p,q)dq = and JW{p,q)dp = How- 

ever, it cannot be identified with a probability distribu- 
tion since it can take negative values. The Wigner func- 
tion has been measured experimentally in atomic sys- 
tems, and such negative values have been reported 25]. 

Although the Wigner function can take negative val- 
ues, it can be shown that coarse graining this function 
over cells of size h always leads to nonnegative values. 
Therefore a smoothing of Q by appropriate functions 
will lead to a function of p and q with no negative values. 
An example of such a function is given by the Husimi dis- 
tribution (see e.g. which uses a Gaussian smoothing. 
A further example using another smoothing function was 
discussed in 14]. 

In the following sections, we will study the evalua- 
tion of such quantum phase space distributions of wave 
functions on a quantum computer. This will be per- 
formed using a specific example, namely the kicked ro- 
tator model. This system corresponds to the quan- 
tization of the Chirikov standard map jl^, 113 ^ — 
n-\- k sin 0] = -\- Tn where (n, 0) are the conjugated 
(action- angle) variables. 

The classical standard map depends only on the pa- 
rameter K = kT. The system undergoes a transition 
from integrability {K = 0) to more and more devel- 
oped chaos when K increases, following the Kolmogorov- 




FIG. 1: Classical phase space distribution for the standard 
map with K = 0.5 (top left), K = 0.9 (top right), K = 1.5 
(bottom left), K — 2 (bottom right). Black is zero probabil- 
ity, white is maximal probability. As initial state we chose a 
uniform distribution on the set — tt < p < — 3/47r, < x < 27r, 
and 1000 iterations of the standard map were performed. 



Arnold-Moser theorem. Chaotic zones get larger and 
larger until the value K = Kg ^ 0.9716... is reached, 
where global chaos sets in, but a complex hierarchical 
structure of integrable islands surrounded by chaotic lay- 
ers is still present. For K ^ Kg^ the chaotic part covers 
most of the phase space. This system has been used for 
example as a model of particle confinement in magnetic 
traps, beam dynamics in accelerators or comet trajecto- 
ries j^. Its phase space is a cylinder (periodicity in ^), 
and since the map is periodic in n with period 27r/T, 
phase space structures repeat themselves in the n direc- 
tion on each cell of size 2it/T. FigQshows one such phase 
space cell for various values of the parameter show- 
ing the different regimes from quasi-integrability (many 
invariant curves preventing transport in the momentum 
direction) to a mixed regime with a large chaotic domain. 

The quantum version of the standard map 26] gives a 
unitary operator acting on the wave function through: 



-ik cos e -iTn^ / 2 



(2) 



where h = -id/dO, h=l, and V^((9 + 27r) = V^(i9). 

The quantum dynamics Q depends on the two param- 
eters k and T, T playing the role of an effective h. The 
classical limit is /c ^ oo, T ^ while keeping K = kT = 
constant. 

This quantum kicked rotator Q is described by quite 
simple equations, making it practical for numerical sim- 
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ulations and quantum computing. Nevertheless, it dis- 
plays a wealth of different behaviors depending on the 
values of the parameters. Indeed, classical dynamics un- 
dergoes a transition from integrability to fully developed 
chaos with intermediate mixed phases between these two 
regimes. Wave functions show complex structures re- 
lated to the classical phase space corresponding to these 
different cases. In addition, for large K where classical 
dynamics is strongly chaotic, quantum interference can 
lead to exponential localization of wave functions. This 
phenomenon is related to the Anderson localization of 
electrons in solids, and therefore enables to study this 
important solid state problem, which is still the subject 
of active research. The kicked rotator can also model 
the microwave ionization of Rydberg atoms ji^ . an d has 
been experimentally realized with cold atoms 29]. For 
all these reasons, it has been the subject of many stud- 
ies, and can be considered as a paradigmatic model of 
quantum chaos. 

In 0, it was shown that evolving a TV-dimensional 
wave function through the map Q can be done with only 
0{\ogN) qubits and 0{{\ogN)'^) operations on a quan- 
tum computer (compare with O(A^logA^) operations for 
the same simulation on a classical computer). Another 
quantum algorithm developed in 9] enables to perform 
the same quantum evolution (albeit approximately) with 
0((logA/')^) operations. This system can therefore be 
simulated efficiently on a quantum computer, and can be 
used as a good test ground for assessing the complexity 
of various quantum algorithms for quantum phase space 
distributions. 

In the following sections, we will study the efficiency of 
various quantum algorithms to obtain various informa- 
tion about the quantum phase space distribution func- 
tions. The simulation of a quantum system on a quan- 
tum computer based on qubits implies that the system 
is effectively discrete and finite. We therefore close the 
phase space in the momentum direction through peri- 
odic boundary conditions. We will concentrate on the 
regime where T = 27r/A^, N being the Hilbert space di- 
mension. This implies that the phase space contains only 
one classical cell, and increasing the number of qubits at 
K constant decreases the effective h keeping the classical 
dynamics constant. Different K values enable to probe 
various dynamical regimes, from integrability to chaos. 
The localization length in this regime becomes quickly 
larger than the system size for small number of qubits, 
thus allowing to explore the complexity of a chaotic wave 
function. Indeed, in the localized regime, the most im- 
portant information resides not so much in such distribu- 
tions, but in the localization properties, and their mea- 
surement on a quantum computer was already analyzed 
in [illl^. 

For such a quantum system on a dimensional Hilbert 
space, the general formalism of Wigner functions should 
be adapted. In particular, it is known that it should 




FIG. 2: Wigner function for the quantum kicked rotator 
with parameters of FigQ = 0.5 (top left), K = 0.9 (top 
right), K = 1.5 (bottom left), K = 2 (bottom right). Here 
T 27r/A, where A = 2''^ with riq 7. The whole Wigner 
function (on a 2 A x 2 A lattice) is plotted. White marks 
positive maximal values, black negative values. Initial state 
is uniformly spread on the set < n < A/8 (corresponding 
to the initial classical distribution in Fig.l) (this state can be 
built efficiently from |n = 0) by rig — 3 Hadamard gates), and 
Wigner function is computed after 1000 iterations of 



be constructed on 2A^ x 2A^ points (see e.g. |3l|). For 
the kicked rotator, the formula for the discrete Wigner 
function is: 



N-l 



w{e,n)=J2 



'-n{m-e/2) 



2N 



-tp{e-myi^{m), (3) 



with e = ^. 



The Wigner function provides a pictorial representa- 
tion of a wave function which can be compared with the 
classical phase space distribution, (see example in Fig|2j, 
although quantum oscillations are present. 



III. MEASURING THE WIGNER 
DISTRIBUTION 

In the first quantum algorithm was set up which 
enables to measure the value of the Wigner function at 
a given phase space point. The algorithm adds one an- 
cilla qubit to the system and proceeds by applying one 
Hadamard gate to the ancilla qubit, then a certain op- 
erator [/(6,n) is applied to the system controlled by 
the value of the ancilla qubit. After a last Hadamard 
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gate is applied to the ancilla, its expectation value is 
< cr^ >= Re[Tr{U{e,n)p)] = 2NW{G,n) where p is 
the density matrix and N = 2"^-? is the dimension of the 
Hilbert space. One iteration of this process requires only 
a logarithmic number of gates. Nevertheless, the total 
complexity of the algorithm may be much larger, since 
measuring < cr^ > may require a very large number of 
measurements. This can be probed only through care- 
ful estimation of the asymptotic behavior of individual 
values of the Wigner function. 

A drawback of the approach of 15] is that it does not 
allow easily further treatment on the Wigner function 
which may improve the total complexity of the algorithm. 
To this aim, the simplest way is to build explicitly the 
Wigner transform of the wave function as amplitudes of 
a register. This enables to use additional tools (ampli- 
tude amplification, wavelet transforms) which may in- 
crease the speedup over classical computation, as we will 
see. 

Such an explicit construction of the Wigner function 
directly on the registers of the quantum computer is in- 
deed possible in the following way. To get the Wigner 
function of U*\ipo) {t iterations of an original wave func- 
tion |?/^o) through we start from an initial state (for 
example in n representation) \ipo)'S)\ipo) = ^fjo^ 

T.^Jo^(^jM = Ez=o^Ej^"o^«iS*l^^)l^^)- ™^ ^^^^^ 
2nq qubits to hold the values of the wave function on 

a A^-dimensional Hilbert space, where N = 2^"?. Then 
we apply the algorithm implementing the kicked rotator 
evolution operator U developed in Q to each subsystem 
independently. This operator can be described as multi- 
plication by phases followed by a quantum Fourier trans- 
form (QFT). The multiplication by phases of each coef- 
ficient keeps the factorized structure. The QFT mixes 
only states with the same value of the other register at- 
tached, and therefore also keeps the factorized form. Let 
us see how it works for one iteration: 

AT-lAT-l N-lN-1 
i=0 j=0 i=0 j=0 

(multiplication by e~*^^^^/^) 

N-l N-1 N-1 N-1 

j=0 i=0 j=0 i=0 

(QFT with respect to rii) 

N-1 N-1 

(multiplication by e~*^^°®^^) 

N-1 N-1 
j=0 i=0 



(QFT with respect to Oi) 

N-1 N-1 

We can thus get U*\ipo)^U*'^\ipQ) by applying the process 
several times. This can be done in a number of gates 
polynomial in Uq (O(tn^) if we use the algorithm of |^ 
for implementing U). 

From such a state it is possible to build efficiently the 
state nW{0^n)\0)\n). Indeed, building the Wigner 
transform can be done through a partial Fourier trans- 
form. To see this, we start from the state in repre- 
sentation, i.e. (g) |V^*)= T.e,e^^(^)^^i^')\^W)' Then 
we add an extra qubit to the first register (needed to get 
values of ^ + between and 2N — 1) and realize the 
transformation * 

Eo,o'mr{o')\ew) ^ Eo,o'mr{o')\e + ew) 

(addition) 

Let us call O = -\- 0'; then the state can be 
written Ee,^' - 6>OV^*((9O|0)|6>O Then we real- 
ize a QFT of the second register only. The re- 
suit is: EeEnlE^'e-'^^'Vl© - 0')rmm\n) 
=2V^J]@ W(e,n)e-^^®/2|e)|n) where 9 varies 
from to 2Ar - 1 and n from to AT - 1. To 
get the Wigner function on a 2A^ x 27V grid, we 
need first to add an extra qubit in the state |0) 
and apply a Hadamard gate to it. If we inter- 
pret it as the most significant digit of n, the result- 
ing state is V2NY:^QY^^~^W{e,n)e-^''^/^\e)\n) + 
^EeEn=^'^(0,^)e"'^^""'^^®/'|e)|n). The fi- 
nal step consists in multiplying by the phases e~~w^^l'^ 
and e~"^^"^~^^®/^, which can be made by application 
of two-qubit gates (controlled phase-shifts). The final 
state is 

2N-1 2N-1 

\^f) = V2Nj2 E W{e,n)\e)\n) (4) 

One can check that the normalization is correct since 
it is known in general that Ee^o^ En=o ^ ^(®' 
1/2N. 

The advantage of this procedure in comparison to the 
one in 15] resides in the fact that individual values of 
the Wigner function are now encoded in the compo- 
nents of the wave function. This is in general a natu- 
ral way to encode an image on a wave function: each 
basis vector corresponding to a position in phase space 
is associated with a coefficient giving the amplitude at 
this location. This way of encoding the Wigner func- 
tion enables to perform some further operations to ex- 
tract information efficiently through quantum measure- 
ments. We will envision three different strategies: direct 
measurements of each qubit, amplitude amplification and 
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wavelet transform. The data of Fig l3l6l will enable us to 
compare these different strategies for different physical 
regimes of the kicked rotator model, with various levels 
of chaoticity. The quantity plotted is the inverse partici- 
pation ratio (IPR). For a wave function = Yl!i=i V^^K)? 
where \i) is some basis, the inverse participation ratio is 
measures the number of significant 
components in the basis The Wigner function veri- 
fies the sum rules = 1 and Wf = . Following 
we are lead by analogy to define the inverse par- 
ticipation ratio for the Wigner function, by the formula 
^ = ^/ {N'^^W^). If the Wigner function is composed 
of N components of equal weights then ^ = 
whereas N'^ components of equal weights (in absolute 
value) l/7V^/2 give ^ = Thus the IPR ^ gives an 
estimate of the number of the main components of the 
Wigner function. 
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FIG. 3: Main plot: scaling of the IPR ^ vs. Uq for the Wigner 
function (empty squares) and for the wavelet transform of 
the Wigner function (full squares). The full straight lines 
represent the law iV^, N ^ 2"^"^ . Here K = 0.5. In the 
inset, the ratio R between IPR of Wigner function and wavelet 
transformed Wigner function is plotted. Parameters, number 
of iterations and initial state are the same as in Fig|21 

To compare classical and quantum computation of this 
problem, we first should assess the complexity of obtain- 
ing the Wigner function on a classical computer. For a 
A/'-dimensional wave function, iterating t times the map 
costs 0{tN log TV) operations. Then getting all values 
of W needs to perform N Fourier transforms, requiring 
0{N'^ log N) operations. The same is true for obtaining 
the largest values of W, if one does not know where they 
are: only the computation of all of them and subsequent 
sorting can provide them. Thus in both cases classical 
complexity is of the order 0{N^ \ogN). This asymptotic 
law changes if one is interested in a single value of the 
Wigner function at some predetermined (9,n) value. In 
this case, only one Fourier transform is actually needed, 
so the classical complexity becomes of order 0{N log N). 

As concerns the quantum computer, we have to clarify 
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FIG. 4: Main plot: scaling of the IPR ^ vs. riq for the 
Wigner function (empty squares) and for the wavelet trans- 
form of the Wigner function (full squares). The full straight 
line represents the law N'^'^^ , while the dashed line repre- 
sents A^^, iV = 2"''^. Here K = 0.9. In the inset, the ratio 
R between IPR of Wigner function and wavelet transformed 
Wigner function is plotted. The full line represents the scal- 
ing A/"*^"^^. Parameters, number of iterations and initial state 
are the same as in Fig|2] 
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FIG. 5: Main plot: scaling of the IPR ^ vs. riq for the 
Wigner function (empty squares) and for the wavelet trans- 
form of the Wigner function (full squares). The full straight 
line represents the law A/"^"^, while the dashed line represents 
AT = 2^^^. Here K = 1.5. In the inset, the ratio R 
between IPR of Wigner function and wavelet transformed 
Wigner function is plotted. The full line represents the scal- 
ing N^'^ . Parameters, number of iterations and initial state 
are the same as in Fig|21 



the measurement protocol to assess the complexity of the 
algorithm. The most obvious strategy consists in measur- 
ing all the qubits after explicit construction of the wave 
function and accumulating statistics until a good pre- 
cision is attained on all values of the Wigner function. 
From Fig l3l6l (empty squares), we can see that in the 
four physical regimes considered, the IPR scales approx- 
imately as . This implies that the Wigner function 
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FIG. 6: Main plot: scaling of the IPR ^ vs. Uq for the Wigner 
function (empty squares) and for the wavelet transform of 
the Wigner function (full squares). The full straight line rep- 
resents the law N^'^ ^ while the dashed line represents N^'^ ^ 
iV = 2''«. Here K = 2. In the inset, the ratio R between IPR 
of Wigner function and wavelet transformed Wigner function 
is plotted. The full line represents the scaling A/"*^'^^, while the 
dashed line represents iV°"^^. Parameters, number of itera- 
tions and initial state are the same as in Fig|21 



is spread out on the components, each term having 
comparable amplitude Wi ~ N~'^l^ . This needs mea- 
surements to get a good precision. The number of quan- 
tum operations is therefore 0{tN'^) {N'^ repetitions of t 
iterations) up to logarithmic factors. This should be com- 
pared with the classical complexity of obtaining all values 
of the Wigner function, or only the largest ones, which 
both are of order 0{N'^ log N). This makes the quan- 
tum method no better than the classical one, albeit the 
quantum computer needs a logarithmic number of qubits 
whereas the classical computer needs exponentially more 
bits (A^ bits versus log TV qubits). This can translate in 
an improvement in effective computational time by for 
example distributing the computation over subsystems 
of qubits, and making simultaneous measurements, but 
this is obviously quite cumbersome. 

Still, it can be remarked that for the values of K for 
which the system is most chaotic, the IPR scales with 
a slightly lower power with a ~ 1.8 — 1.9. If this 
is verified asymptotically, then the quantum algorithm 
need only 0{tN^) operations, and a small gain of N'^~^ 
is realized. It is interesting to note that if the classi- 
cal algorithm to compute the evolution of the map were 
more complex, i.e. of order 0{tN^) with /c > 2, then 
in this case the quantum algorithm will be better by an 
additional factor of N^~'^. 

The phase space tomography method of requires 
to measure < > of an ancilla qubit, with < cr^ >= 
A^VI^(6,n). Thus < >- 7V-^/^ a value which re- 
quires measurements to be reasonably assessed. This 
should be compared with the classical cost of obtaining 



the value of the Wigner function at a predetermined lo- 
cation, which is of order 0(A^log7V). Again, the method 
is not better than the classical one, although it uses a 
logarithmic number of qubits which can translate into 
an improvement in effective computational time by dis- 
tributing the quantum computation. Similarly, when the 
IPR scales as with a < 2, then the quantum algo- 
rithm is better by a factor of N'^~^. For other maps for 
which classical simulation is of order 0(tN^) with A: > 2, 
the quantum algorithm will be better by an additional 
factor of N^~'^. 

It is possible to use coarse-grained measurements in or- 
der to decrease the number of measurements of the wave 
function To this aim, one measures only the first n/ 
qubits with n/ < Uq {N = 2"^^). This gives the integrated 
probability inside the 2^"^^ cells (sum of 2^^-?"^"^/ proba- 
bilities |VK(0, n)p)) in a number of measurements which 
scales with the number of cells and not any more with the 
number of qubits. This is possible if the wave function of 
the computer encodes the full Wigner function in its com- 
ponents, as in the algorithm exposed above. In principle, 
the complexity is 0{2'^^f) and a gain compared to classi- 
cal computation can be obtained. There is a possibility 
of exponential gain with this strategy, since by fixing n / 
and letting riq increase, measuring the integrated prob- 
ability becomes polynomial in Uq. Still, the precision 
is also polynomial, and it is possible that semiclassical 
methods enable to get such approximate quantities since 
with Uq ^ oo the value of h becomes smaller and smaller 
and the system is well approximated by semiclassical cal- 
culations. If this holds, the advantage of quantum com- 
putation may be less spectacular. 

A similar method can be applied to the phase space 
tomography method of p^, but with a different re- 
sult. In 16] it is explained that one can compute av- 
erages of Wigner function on a given rectangular area 
by using an ancilla qubit. The process gives < >= 
2N^W{Q^n)/Np, where Np is the number of points 
over which the summation is done. Note that contrary 
to the previous discussion, the sum is over W and not 
In this case, the normalization constant Np makes 
the method comparable to direct phase space tomogra- 
phy of one value of the Wigner function at one phase 
space point. With this technique, there is no additional 
gain in adding up components. 

A more refined strategy uses amplitude amplification 
It is a generalization of Grover's algorithm 3]. The 
latter starts from an equal superposition of N states, 
and in operations brings the amplitude along one 
direction close to one. Amplitude amplification increases 
the amplitude of a whole subspace. If P is a projector 
on this subspace, and V is the operator taking |0) to 
a state having some projection on the desired subspace, 
repeated iterations ofV{I-2\0){0\)V-\l-2P) on ^10) 
will increase the projection. Indeed, if one write V\0) = 
PF|0) + (/ - P)F|0), the result of one iteration is to 
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rotate the state toward PV\0) staying in the subspace 
spanned by PV\0) and (/ - P)V\0). If a = \PV\0)\^, 
one can check that after one iteration the state is (4a^ — 
3)PV'|0) + (4a2 - 1)(/-P)V'|0), with a component along 
{I -P)V\0) decreased by 4a^ 

If V is chosen to be UwignerU^ (where Uwigner builds 
the Wigner transform), and P to be a projector on the 
space corresponding to a square of size Nd x Njj^ the 
process of amplitude amplification will increase the total 
probability in the square, keeping the relative amplitude 
inside the square. This acts like a "microscope" , increas- 
ing the total probability of one part of the Wigner func- 
tion but keeping the relative details correct. The total 
probability in a square of size Nd x Nd^ following the 
results shown in Fig. 3-6, should be of the order N'^/N'^. 
Amplitude amplification will therefore need N/Njj iter- 
ations to bring the probability inside the square close to 
one. Then according to Fig 3-6 measurements are 
needed to get all relative amplitudes with good preci- 
sion. Total number of quantum operations is therefore 
0{tNDN) (up to logarithmic factors). This should be 
compared to the number of classical operations, 0{tN) 
for the evolution of the wave function, and 0{NdN) 
for computing the Wigner function (construction of the 
Wigner function on a square of size N^^ needs only Nd 
Fourier transforms of N dimensional vectors). Both com- 
putations are therefore comparable for low K. When the 
scaling N'^ {a < 2) for the IPR of the Wigner function 
is verified, then measurements are enough to get the 
Wigner function on a quantum computer, and a small 
gain of A^^~" is present for the quantum algorithm. If 
the classical algorithm to compute the evolution of the 
map were more complex, i.e. of order 0{tN^) with k > 1 
(note the difference with the previous case where k > 2 
was needed) , then in this case the quantum algorithm 
will be better by an additional factor N^~^ if Nd and t 
are kept fixed. 

Our last strategy uses the wavelet transform. This 
transform 0, is based on the wavelet basis, which 
differs from the usual Fourier basis by the fact that each 
basis vector is localized in position as well as momen- 
tum, with different scales. The basis vectors are ob- 
tained by translations and dilations of an original func- 
tion and their properties enable to probe the different 
scales of the data as well as localized features, both in 
space and frequency. Wavelet transforms are used ubiq- 
uitously on classical computers for data treatment. Al- 
gorithms for implementing such transforms on quantum 
computers were developed in j^O, 0, US , and were 
shown to be efficient, requiring polynomial resources to 
treat an exponentially large vector. Effects of imperfec- 
tions on a dynamical system based on the wavelet trans- 
form were investigated in |2^. In the present paper, we 
implemented the 4-coefficient Daubechies wavelet trans- 
form (P)^^^), the most commonly used in applications, 
and applied it to the two-dimensional Wigner function 



The results in Fig I3I6I show that the IPR of the wavelet 
transform of the Wigner function scales as ^ with P de- 
creasing from /3 ~ 2 to ^ 1.4 when the chaos parameter 
K is increased (for K = 0.5, with low level of chaos, the 
wavelet transforms yield a compression factor of order 10, 
but no visible asymptotic gain) . This means that getting 
the most important coefficients in the wavelet basis needs 
only measurements. The quantum algorithm for get- 
ting them needs only 0{tN^) operations. On a classical 
machine, the slowest part is still the computation of the 
Wigner function, which scales as 0{N'^). Therefore at 
fixed t the gain is polynomial, of order 0{N'^~^). How- 
ever, recovering the coefficients of the original Wigner 
function needs to use a classical wavelet transform which 
needs 0{N'^) operations. Still, the wavelet coefficients 
give information about the hierarchical structures in the 
wave function, so obtaining them with a better efficiency 
gives some physical information about the system. 

Therefore, as concerns the quantum computation of 
the Wigner function, it seems a modest polynomial gain 
can be obtained by different methods, especially in the 
parameter regime where the system is chaotic, the most 
efficient method being the measurement of the wavelet 
transform of the distribution, although the interpretation 
of the results is less transparent. 



IV. MEASURING HUSIMI FUNCTIONS 

As already noted in Section II, the Wigner function is 
comparable to a classical phase space distribution, but 
can take negative values. It is known that it becomes 
non-negative when coarse-grained over cells of size h. 
One way to do this coarse-graining is to perform a con- 
volution of the Wigner function with a Gaussian, giving 
the Husimi distribution 11]: 

pH{Oo,no) = K<^(^o,no)lV^)l^ (5) 

where (^, n) = A^^e'^^-^^^" /^-"-''^^\n) is a 

Gaussian coherent state centered on (^o, ^o) with width a 
{A is a normalization constant). An interesting quantum 
algorithm was proposed in [iQ to compute this distribu- 
tion, based on phase space tomography. It uses a rel- 
atively complicated subroutine which builds an approx- 
imation of coherent states on a separate register. This 
method is similar to the Wigner function computation 
through an ancilla qubit analyzed in the preceding sec- 
tion, and gives comparable results. 

In 14], a very fast quantum algorithm was proposed 
to build a modified Husimi function, which is defined by: 

P^S\(^o,m) = \{<p^^l^J^)f (6) 
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where 4>tLo^(^^ = iVN^") E^ltr"' ^-^''""In) is a 
modified coherent state centered on {Oq^tiq). The convo- 
lution is not made any more with a Gaussian function, 
but with a box function of size in momentum. This 
imphes a very good locahzation in momentum, but in 
contrast in the angle representation the amplitude de- 
creases only as a power law since the Fourier transform 
of the box function is the function ^Hi^. 

X 

This transform can be evaluated quite efficiently on a 
quantum computer without computing the Wigner func- 
tion itself. Indeed, as shown in 14], it can be computed 
by applying a QFT to the ffist half of the qubits. This 
partial Fourier transform uses + 1) quantum el- 

ementary operations to build from a wave function |?/^) 
with N = 2^^^ components the state 

|^j^) = ^F(^,n)|0)|n) (7) 

where and n take only a/TV values each and \H{0^ n)p is 
the modified Husimi function Performing the same 
task on a classical computer needs 0(A^ log(A^)^) opera- 
tions. We will concentrate on this method to compute 
Husimi functions, since it seems to be the most sim- 
ple and easy to implement, and gives a good picture of 
the wave function as can be seen in the implementations 
made in ii,l3i|. 

In Fig[71 we show the result of performing the evo- 
lution on a wave packet in A/'-dimensional Hilbert 
space for four different values of and then applying 
the partial Fourier transform. The result is an array of 
X points, each point representing an average over 
~ N neighboring values of the Wigner function. The fig- 
ure shows that this transformation allows to obtain effi- 
ciently a positive phase space distribution which can be 
compared with the classical distributions, for example in 
FiglH 

In Fig l8lll( we show the IPR of the result of this trans- 
form and of an additional wavelet transform of this func- 
tion. The data show that with this modified Husimi dis- 
tribution the compression of information is much better 
than in the case of the Wigner function. 

Indeed, in all four parameter regimes considered, the 
IPR of the function scales as ^ with 0.5 < 7 < 0.7. 
This means that the most important components of the 
modified Husimi distribution can be measured with ~ 
quantum measurements. Thus on a quantum computer 
the whole process of evolving the wave function up to 
time t, transforming it into the modified Husimi distri- 
bution and measuring it needs 0{tN^) operations. On 
the contrary, a classical computer will need 0{tN) opera- 
tions for the evolution, and 0{N) for the modified Husimi 
transform (up to logarithmic factors). Thus for the sys- 
tem 0, computation of the modified Husimi transform 
is more efficient on a quantum computer (including mea- 
surement) than on a classical one. This gain would dis- 




FIG. 7: (Color online) Modified Husimi function for the 
quantum kicked rotator with K = 0.5 (top left), K = 0.9 
(top right), K = 1.5 (bottom left), K = 2 (bottom right). 
Here T = 27i/N, where N ^ 2""^ , with Uq = 16. The function 
is plotted on a lattice of >/]V x a/]V and each point is the 
average of N points. Initial state is the same as in Fig|21 
(corresponding to the initial classical distribution in Fig^, 
and the function is computed after 1000 iterations of Red 
(gray) is maximal value, blue (black) minimal value. 



appear if the transform had IPR ~ TV. 

As in the preceding section, one can use coarse-grained 
measurements in order to increase the probability. Again, 
this gives the integrated probability inside the cells in a 
number of measurements which scales with the number 
of cells, with the same drawbacks than in section HI. 

If we use amplitude amplification, the gain is even 
better. Amplitude amplification will need ^/N/Nd it- 
erations to bring the probability inside a square of size 
^/N]j X y/N]j close to one. Then according to Fig. 8- 
11 measurements are needed, with 0.5 < 7 < 0.7. 
The total number of quantum operations is therefore 
0{tVNN]^~^^'^). Classically, we stih need 0{tN) opera- 
tions for the evolution, and O^-s/Ny/ Nd) for the trans- 
form (up to logarithmic factors). Thus for small Nd 
a quadratic gain is achieved. Interestingly enough, this 
gain persists in the case where the IPR of the modi- 
fied Husimi function is ~ A/", even though the previous 
method then will not give any gain. Since the IPR cannot 
be larger than A/", this means that with amplitude am- 
plification the quantum computer is in general at least 
quadratically faster at evaluating part of the modified 
Husimi function than any classical device. If the classi- 
cal algorithm to compute the evolution of the map were 
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FIG. 8: Scaling for the IPR ^ vs Uq for the function H{0, n) in 
Q, with parameters K = 0.5 and T = 27r/iV, iV = 2'''^ . Main 
plot: empty squares represent the IPR of H(0,n) function , 
the full squares represent the IPR of the wavelet transform of 
the modulus of H{0,n). The fuh line is iV°-^. In the inset, 
the ratio R between IPR of H(0,n) and wavelet transform 
of \H{0,n)\ is plotted for different rig, the fuh line is N^'^. 
Number of iterations and initial state are the same as in Fig|71 
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FIG. 10: Scaling for the IPR i vs for the function i7 (6*, n) 
in O, with parameters K 1.5 and T 27r/A^, AT = 2''^. 
Main plot: empty squares represent the IPR of H(0,n), the 
full squares represent the IPR of the wavelet transform of the 
modulus of H{0,n). The full line is N^'^, the dashed line 
is N'^'^ . In the inset, the ratio R between IPR of H(0,n) 
and wavelet transform of \H{0,n)\ is plotted for different Uq. 
Number of iterations and initial state are the same as in FigQ 
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FIG. 9: Scaling for the IPR ^ vs riq for the function H{6,n) 
in 0, with parameters K 0.9 and T 27t/N, N = 2""^ . 
Main plot: empty squares represent the IPR of H(0,n), the 
full squares represent the IPR of the wavelet transform of the 
modulus of H(6,n). The full line is N^'^, the dashed line is 
N^-^. In the inset, the ratio R between IPR of H{6,n) and 
wavelet transform of \H{0,n)\ is plotted for different n^, the 
full line is N^'"^ . Number of iterations and initial state are the 
same as in FigQ 




FIG. 11: Scaling for the IPR ^ vs Uq for the function H{0, n) 
in Q, with parameters = 2 and T = 27r/Ar, = 2""^ . 
Main plot: empty squares represent the IPR of H(0,n), the 
full squares represent the IPR of the wavelet transform of the 
modulus of H{0,n). The dashed line is N^'^ . In the inset, 
the ratio R between IPR of H(0,n) and wavelet transform 
of \H{0,n)\ is plotted for different rig, the full line is N^'^ . 
Number of iterations and initial state are the same as in FigjTl 



more complex, i.e. of order 0{tN^) with /c > 1, then in 
this case the quantum algorithm will be better by a fac- 
tor A^^~^/^ if Nd and t are kept fixed, making the gain 
even larger, but still polynomial. 

We also analyzed the use of the wavelet transform to 
compress these data and minimize the number of mea- 
surements. At this point a slight complication appears. 
In the previous section, individual amplitudes of the wave 
function in were actual values of the Wigner function. 



so performing a quantum wavelet transform of was 
equivalent to a wavelet transform of the Wigner function. 
In the case at hand, the wave function of the quantum 
computer is such that the modulus square of its compo- 
nents give the modified Husimi distribution. One can 
perform a quantum wavelet transform of this wave func- 
tion, with real and imaginary parts for all coefficients, 
which gives the wavelet coefficients of a complex- valued 
distribution whose square is the modified Husimi distri- 
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bution. It is not clear how to interpret the resulting 
coefficients, and anyway our data have shown that this 
process does not decrease the IPR (data not shown), thus 
making it an inefficient way of treating such data. How- 
ever, Figures 8-11 show that if one take the modulus of 
the wave function, then the IPR of the wavelet trans- 
form of this function is quite small, scaling as 0{N^)^ 
with S ^ — 0.2. So the modified Husimi function itself 
is well compressed by the wavelet transform. It is the 
phase of H{0^ n) in (0) which, although irrelevant for the 
Husimi functions, prevents compression by the wavelet 
transform. To use efficiently the wavelet transform, we 
therefore need to get rid of the phases, i.e. construct a 
wave function whose components are the moduli or mod- 
uli square of the preceding wave functions. 

Such a wave function can be prepared by starting from 
two initial wave packets on two separate registers ^ 
l-^J), and as in the previous section make them evolve 
independently to get U^\^lJ{)) Then a partial 

Fourier transform is applied independently to both reg- 
isters, yielding Y.H{0,n)H{0' ,n'Y\e)\e')\n)\n') . Then 
amplitude amplification should be used to select the di- 
agonal components, yielding ^\H{0.,n)\^\9)\n) . These 
components represented a probability N/N'^ = 1/N of 
the full original wave function, thus this process costs 
0{t^/N) operations up to logarithmic factors. This pro- 
cedure gives us a final wave function whose components 
are now the modified Husimi function itself, without 
the irrelevant phases. We can now apply the quantum 
wavelet transform to this wave function. Afterward, mea- 
suring the main components of the wave function should 
need only 0{N^) quantum measurements. The cost of 
the total procedure is therefore 0{tN^~^^^'^) quantum op- 
erations, whereas classical computation will cost 0{tN) 
operations. Obtaining the main wavelet components of 
this modified Husimi distribution is therefore more effi- 
cient on a quantum computer than on a classical one, 
albeit the gain is still polynomial. 

V. STANDARD IMAGES 

The investigations in the previous sections show that 
computation of quantum phase space distributions can 
be more efficient on a quantum computer than on a clas- 
sical device. An usually polynomial gain can be obtained 
for the whole process of producing the distribution and 
measuring its values. These phase space distributions 
are in effect examples of two-dimensional images. It is 
interesting to explore these questions of efficiency of im- 
age processing on a quantum computer in a more general 
setting. 

FigEl shows four examples of classical images which 
we can use as benchmark to test different strategies of 
processing them. The top left image (the girl) is a stan- 
dard example used in the field of classical image process- 




FIG. 12: Images analyzed in this section. Top: girl image 
(left) and New York City picture (right). Bottom: galaxy 
image taken from NASA website (left) and a fractal picture 
built on the purpose of studying image compression (right). 



ing. Top right is a aerial view of New York City, bot- 
tom left an astronomical photograph, and bottom right 
an artificially built picture with fractal-like structures. 
They represent diverse types of images that can be pro- 
duced and processed for various purposes. We will sup- 
pose in the following that these black and white pictures 
are encoded on a quantum wave function in the form 
^ — y ^xy\x)\y) where x^y are indexes of N'^ pix- 
els and axy are the amplitudes on each pixel (positive 
number). Of course, contrary to the previous examples, 
we do not know how to produce in an efficient way such 
types of wave functions. We therefore concentrate on the 
problem of extracting information efficiently from such a 
wave function once it has been produced. 

Fig^l permits to analyze two of the strategies pre- 
cedingly developed. The IPR of the different images are 
shown to scale like A^^, implying that direct measure- 
ment of all qubits will need 0{N'^) measurements to get 
the most important components (since these components 
scale also like 0{N'^)). As in the case of the Wigner 
function, coarse-grained measurements are possible, and 
require a number of measurements proportional to the 
number of cells. This is more efficient, at the price of 
losing information on scales smaller than the cell size. 

The use of amplitude amplification on a small part of 
the picture (polynomial in Uq) enables to bring this part 
to a probability close to 1 in 0{N) Grover-like iteration. 
So if one is interested in details of the picture at a specific 
place predetermined, this strategy is more efficient than 
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FIG. 13: (Color online) Scaling for the IPR ^ vs Uq for the im- 
ages with different resolutions. Full symbols: ^ after wavelet 
transform. Empty symbols: ^ for the original images, with 
different resolution (from 32 x 32 to 2048 x 2048). Squares 
refer to the girl image, circles to the New York City image, 
triangles to the galaxy image, and diamonds to the fractal 
image. The dashed hne is the law iV^ with = 2""^. The 
original images are 8-bit gray scale images. They are encoded 
in the wave function from which the IPR is computed. 



the direct measurement. Of course, precise efficiency of 
the quantum process compared to classical methods will 
depend on the relative complexity of the classical and 
quantum image production, which probably varies with 
the problem. 

The full symbols in FigElgive the IPR of the wavelet 
transform of the image. That is, the image is encoded in 
a quantum wave function as previously, and a quantum 
wavelet transform is applied. The resulting wave func- 
tion displays an IPR which grows slowly with Uq. Actu- 
ally, data from FigEl are compatible with a polynomial 
growth with Uq of the IPR. This would indicate that the 
wavelet transform is very efficient in compressing infor- 
mation from standard images. Obtaining the main com- 
ponents of the wavelet transform would demand polyno- 
mial number of measurements compared to an exponen- 
tial one for the original image wave function. This can 
transfer to an exponential gain in the full process if the 
image can be encoded also in a polynomial number of 
operations in riq. 

In Fig tT^ a different strategy is studied. Namely, in 
analogy with the MPEG standard for image compression, 
the image is decomposed into many tiles, and each tile 
is independently wavelet-transformed. This procedure is 
tested in the case where tiles are of size a/TV x FigtT^ 
shows that although the final IPR grows more quickly 
with Uq than in the case of Fig tT3l the IPR seems asymp- 
totically to be smaller again than with the full image wave 
function. Data from Fig^^are compatible with an IPR 
growing like O(A^), implying that the number of mea- 
surements is the square root of the one for the full wave 



FIG. 14: (Color onhne) Scahng for the IPR ^ vs Uq for the 
images with different resolutions, the tiling method (see text) 
is used, with tiles of size ^ ^/N x >/]V Full symbols: ^ af- 
ter wavelet transform. Empty symbols: ^ for the original im- 
ages, with different resolution (from 32 x 32 to 2048 x 2048). 
Squares refer to the girl image, circles to the New York City 
image, triangles to the galaxy image, and diamonds to the 
fractal image. The dashed line is the law with N = 2^*^, 
the full line is the law N. The original images are 8-bit gray 
scale images. They are encoded in the wave function and the 
IPR is computed from the latter. 
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FIG. 15: Comparison IPR / entropy for the girl image of 
Fig ll2l Full symbols are for IPR, empty symbols for 2*^, where 
S is the entropy. Squares and circles are for the wavelet trans- 
form, diamonds and triangles for the original image. Data for 
the three other images of Fig ll2l give the same result. 



function. This suggests a polynomial speed up with this 
method. We note that a similar strategy for a quantum 
sound treatment was discussed in |3^ . 

FigEl enables to confirm the preceding results which 
use the IPR. Indeed, an alternative quantity to quantify 
the spreading of a wave function on a predetermined ba- 
sis is the entropy. For a A^-dimensional wave function Itp) 
with projections on a basis \(j)j) given by Wj = 
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the entropy is defined by = — Wj log2 Wj . It takes 
values from S = {iIj = (j)j for some j) to 5* = log2 N 
— Both IPR and 2*^ give an estimate 

of the number of components of the wave function. The 
data show that although both quantities are different, 
they show a similar behavior with as do their wavelet 
transform, confirming that the preceding results are ro- 
bust. 




FIG. 16: Image reconstruction from Monte Carlo sampling 
and quantum wavelet sampling. Top: exact girl image (left) 
and sampled with 2500 Monte Carlo points (right). Bottom: 
reconstruction after sampling with 2500 measurements in the 
wavelet basis after full (left) and tiled (right) wavelet trans- 
form. The images are 128 x 128 i.e. 16000 points in total. 



some information about it that can be obtained with a 
small number of measurements. 



VI. CONCLUSION 

In this paper, we have analyzed and numerically tested 
the quantum computation of Wigner and Husimi distri- 
butions for quantum systems. Two methods of computa- 
tion for the Wigner function, one original to this paper, 
were considered. We studied different strategies to ex- 
tract information from the wave function of the quantum 
computer, namely direct measurements, coarse-grained 
measurements, amplitude amplification and measure of 
wavelet- transformed wave function. For the Wigner func- 
tion, the largest (polynomial) gain is obtained through 
the use of the wavelet transform, although other meth- 
ods might yield a smaller gain in the chaotic regime. For 
the Husimi distribution, the gain is much larger, although 
it is still polynomial, and increases with the use of ampli- 
tude amplification and wavelet transforms. At last, the 
study of real images show that the wavelet transform en- 
ables to compress information and therefore to lower the 
number of measurements in the quantum case, although 
a lot of information is lost in the process. 
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the EC RTN contract HPRN-CT-2000-0156 and by the 
project EDIQIP of the IST-FET program of the EC. 



The preceding discussion gives some numerical argu- 
ments suggesting that main components of the wavelet 
transform can be obtained more efficiently than the im- 
age itself. This gives information on the patterns present 
in the picture, and can be considered as an information 
in itself. It is also worth studying how much information 
about this original image is present in these main compo- 
nents of the wavelet transform. FigtT^ shows an attempt 
of reconstruction of one image from these main compo- 
nents only. The results displayed on this figure show 
that although some features are distinguishable with this 
technique (better than with the Monte-Carlo sampling), 
a lot of information from the original figure has been lost. 
It is possible that better results are obtained for larger 
system sizes, but this regime cannot be reached by our 
classical numerical simulations. Still, even if the largest 
wavelet coefficients by themselves are not enough to give 
a good approximation of the original image, they bring 
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